Robust morphogenesis by chaotic dynamics

This research illustrates that complex dynamics of gene products enable the creation of any prescribed cellular differentiation patterns. These complex dynamics can take the form of chaotic, stochastic, or noisy chaotic dynamics. Based on this outcome and previous research, it is established that a generic open chemical reactor can generate an exceptionally large number of different cellular patterns. The mechanism of pattern generation is robust under perturbations and it is based on a combination of Turing’s machines, Turing instability and L. Wolpert’s gradients. These results can help us to explain the formidable adaptive capacities of biochemical systems.

Reaction-diffusion model defined by the system of equations (1) is considered by [1]. To simplify the statement we restrict ourselves to the case of the following two-component system: where u = u(x, y, t) and v(x, y, t) are unknown functions defined on Ω ×{t ≥ 0}, Ω is the strip (−∞, ∞) × (0, 1) ⊂ R 2 , D v and d u are diffusion coefficients for reagents u and v, respectively, η(x, y) and ζ(x, y) are continuous functions, which can be interpreted as external sources independent of u, v. We assume The initial conditions are given by v(x, y, 0) = v 0 (x, y), u(x, y, 0) = u 0 (x, y), where u 0 , v 0 ∈ C 0 (Ω). To simplify the problem, we set the periodic boundary conditions v(x, y, t) = v(x + 2π, y, t), u(x, y, t) = u(x + 2π, y, t), assuming that all the functions ζ, η, u 0 and v 0 are 2π -periodic in x. At the boundaries y = 0 and y = 1 we set the zero Dirichlet conditions for v, i.e., v(x, 1, t) = v(x, 0, t) = 0 (A6) and the zero Neumann conditions for u: 1 Dynamical complexity Following [1] let us consider system (2), (3) with d u << D v . Our results hold under fairly general conditions to f , h, which we formulate below.

Assumptions to f and h
Suppose that (Bi) There exist u * and v * such that and h u (u * , v * ) ̸ = 0; (A2) (Bii) The critical point u * , v * is non-degenerate, i.e., the Hessian H f of the function f at the point (u * , v * ) defined by Let us consider certain examples. In many chemical and biological applications we are dealing with quadratic f and g, for example, Moreover, for chemical or biological applications it is natural to assume that reagent concentrations stay positive if they are positive at the initial moment. This means that the positive cone R 2 > is invariant under the corresponding semiflow. Then c 12 > 0 and, therefore, conditions (A10) reduce to c 12 > 0, f 12 < 0 and c 11 > 0.
In [1] the Brusselator is considered, which gives an example with cubic nonlinearities and which is an important theoretical model of an autocatalytic reaction. This system exhibits Turing instability and describes complicated spatial patterns [2]. It can be shown that our assumptions Bi and Bii hold for this classical model [1].
Note that conditions Bi and Bii exclude monotone and gradient systems, which do not exhibit a stable chaotic large-time behavior. In fact, for gradient systems D v f = g u thus conditions (A1) and (A2) cannot be fulfilled simultaneously.
To see a connection between assumptions Bi, Bii and conditions of monotonicity, let us consider the following class of systems with the parameters λ, d, D, η(·, ·) and ζ(·, ·): under boundary and initial conditions (A4) -(A7). Let consider λ i , d u , D v and functions η, ζ as parameters P of this family of the semiflows. Then the initial boundary value problems (IBVPs) are defined by eqs. (A11), (A12), boundary conditions (A5), (A6), (A7 ) and initial conditions (A4) generate a family of local semiflows S t (P ), where each semiflow depends on the parameters P . Since the coefficients λ i are in the parameter list, we can remove the first condition (A1). Indeed, for any u * , v * we can set λ 1 = f u (u * , v * ). Therefore, we conclude that the family of semiflows defined by IBVPs (A11), (A12), (A4) -(A7) has the property of universal dynamical approximation (see SM and [3]) under the following conditions: there exists a (u * , v * ) such that or there exists a (u * , v * ) such that Note that conditions (A13) and (A14) admit a transparent chemical interpretation: one of the reagents (u, v) is neither an activator nor an inhibitor for the other reagent.

Theorem on attractor complexity
Let us formulate the main result on the attractor complexity. To this end, let us remind Theorem on the Persistence of hyperbolic sets (for definitions of topological orbital equivalency, hyperbolic dynamics and structural stability see [4]).
Roughly speaking, hyperbolic dynamics (which may be periodic or even chaotic), is stable under small perturbations: if Γ is a hyperbolic set for a flow S t on a compact smooth finite dimensional manifold, then for any C 1 -close flow S t there exists a hyperbolic setΓ close to Γ and carrying the same dynamical behavior. The last fact means that there exists a homeomorphism h close to identity, which transforms trajectories of S t on Γ onto trajectories ofS t onΓ (this homeomorphism performs orbital topological equivalency). Since h is close to identity, the distance dist(Γ,Γ ) between Γ andΓ is small. Recall that for two subsets A, B of a Banach space B with a norm || || the distance dist(A, B) is defined by Theorem I Suppose assumptions Bi and Bii are satisfied. Then the family of the local semiflows S t P generates all finite-dimensional hyperbolic dynamics (up to orbital topological equivalencies) as we vary the parameters P of the IBVP defined by eqs. (2), (3) and conditions (A4) -(A7).
The solutions u, v can be described explicitly. Namely, there are functions u 0 (x, y), v 0 (x, y) and θ i (x, y), V i (x, y), i = 1, . . . , N such that where X i (t) are functions of t, γ > 0 is a small parameter, andṽ, u 1 are small corrections such that ||ṽ|| α , ||ũ|| α < cγ s0 , with s 0 > 1 and α ∈ (3/4, 1). Here ||f || α denotes the standard norm in fractional Banach spaces associated with our problem (for a definition of these spaces see, for example, [5]): The admissible choice of the numbers N is defined by on d u and D v : Then, as a result of time rescaling τ = γ s t with a s > 0 and removing small correction terms, we obtain the following system of differential equations with quadratic nonlinearities for X i (t): where a constant c does not depend on small parameters, and a square N × N matrix M with entries M ij = M ij (η, ζ) depends, in a complicated way, on the functions η and ζ. Nonetheless, one can prove the following claim, which is a key point in the proof of Theorem I [1].
Control Lemma I.For each δ > 0 and each square N × N matrix W ij there exist smooth η and ζ such that where [1, N ] denotes the set {1, 2, ..., N }. One can show [6] that when we change the number N and the matrix M, the quadratic systems (A18) realize a dense set in the space of all finite-dimensional smooth vector fields defined on unit balls, i.e., the family of systems (A18) enjoys the property of universal dynamical approximation (UDA) (about UDA, see subsection 1.7, and [3]). This fact implies Theorem I, which asserts that the large time dynamics of X i (t) can be complicated when the problem parameters are adjusted in a special way. Consider systems (A11), (A12) under boundary and initial conditions (A4) -(A7). For these IBVPs, we obtain, as a corollary of Theorem I, that either they generate monotone semiflows, or these systems can generate all structurally stable dynamics when we vary their parameters. In fact, if condition (A13) is not fulfilled, then f v does not change its sign. On the other hand, if (A14) does not hold, g u does not change its sign. Therefore, if conditions (A13) and (A14) are not satisfied, both f v and g u do not change their sign and therefore, IBVP is defined by eqs. (A11), (A12) under conditions (A4) -(A7) generates a monotone semiflow [7]. Finally, we obtain that either our model generates monotone semiflows only, or these systems can generate all structurally stable dynamics when we vary their parameters. Note that Theorem I can be extended on the case of dimension d = DimΩ > 2, however, it is not known whether Theorem I is true for d = 1.

Chemical-biological interpretation
Theorem I can be explained by non-formal arguments. The spatially heterogeneous sources η(x, y) and ζ(x, y) create a non-trivial spatially heterogeneous non-perturbed pattern (u 0 (x, y), v 0 (x, y)), and for all t > 0 the solutions (u(x, y, t), v(x, y, t)) stay in a small neighborhood of that pattern. Inside this neighborhood, the dynamics S t (P ) is weakly nonlinear: it is governed by an evolution equation, which involves a linear operator A and small nonlinear terms [1]. Adjusting appropriate η, ζ we obtain that the spectrum of A has a remarkable structure: it contains N of eigenfunctions ψ i = (θ i , V i ) tr with small eigenvalues and all the remaining spectrum lies in a negative left half-plane and it is separated by a non-small barrier b > 0 from away the imaginary axis. On the physical language, we have N slow modes and an infinite number of fast modes. It is well known that the large-time dynamics of such slow-fast systems is captured by dynamics of slow modes [8]. The crucial fact is that we can control the number of slow modes N and a small interaction between those mode amplitudes X i (t) via the spatially heterogeneous sources η(x, y) and ζ(x, y) (see Control Lemma I formulated above). This control property, which holds only for dimensions d ≥ 2, allows us to prove that our semiflow S t (P ) has the property of universal dynamical approximation.
This picture can be translated into biological language. Consider external sources η(x, y) and ζ(x, y) in (2), (3) as spatial concentrations of certain morphogens. To obtain a complicated attractor, we need complicated functions η(x, y) and ζ(x, y). We can thus consider those functions as carriers of the positional information that permits us to create complicated spatiotem-poral patterns. Moreover, we note that in order to have a complex pattern u(x, y, t), v(x, y, t) we need small coefficients d u and D v >> d u , i.e., one reagent should diffuse much faster than the other one. Under this condition, the positional information can be transformed into a spatiotemporal structure with complex large-time behavior. So, our mechanism to generate complex patterns is based on the Wolpert concept of positional information [9] and the Turing approach [10]. Note that the idea to use spatial heterogeneity was proposed still in [10].

Generation of a number of local attractors
In this subsection, we state a new result with respect to [1]. Consider the following relation for M ij obtained in [1]: where b m = 1 − (−1) m , the coefficientsĝ mj depend on the problem parameters P and we can adjust the values of those coefficients in an arbitrary way as well as the pointsx 1 < ... <x N such that where ϵ > 0 is a small parameter. We set N = N cM , whereM and N c are integers and decompose the pointsx j into N c groups G 1 , ..., G Nc . The points of k-th group G k we denote by x i,k . The key idea in the choice of these groups and the points is that the points inside the same group are close while the groups are separated. Namely, as ϵ → 0. Moreover, we set where ⌊x⌋ stands for the floor of x. Then for all i, j ∈ [1,M ] one has as ϵ → 0. Then we obtain that system (A18) can be decomposed into N c independent (up to exponentially small smooth corrections) systems The matrices M (k) ij depends on the parameters P and for them an analog of control Lemma holds. Therefore each system (A26) realizes a dense subset of the set of all vector fields defined on the balls {||X (k) || < R}, where X (k) is the vector with components X j,k , j = 1, ...,M and R is an appropriate constant (on realizations of vector field see Appendix ). Therefore, we can chooseM and M (k) ij such that k-th system (A26) realizes a structurally stable smooth dynamical system of the radius R 0 > 0 centered at q = 0 and having two local attractors. The first one is a rest point {0} and the second one is a hyperbolic chaotic attractor A k . Note that each sequence A 1 A 2 ...A Nc , where A k is either {0}, or A k corresponds to a local attractor of the system (A18). This shows that the system (A18) has at least 2 Nc local attractors, which are Cartesian products of local attractors for subsystems. The number N c admits estimate N c ≥ M a . This observation implies the assertion (i) of Claim I for two-component systems when m = 2 and thus for all m ≥ 2.
To show assertion ii, we can consider a system composed ofm = ⌊m/2⌋ independent subsystems of the form (1). If each subsystem has M local attractors, then the subsystem has Mm local ones. Let us note that weak interactions between subsystems do not violate this property. This observation finishes the proof of II.

Morphogenesis algorithm
We assume that j-th cell occupies a subdomain Let us describe now how the cell pattern can be produced in our model. Let us consider organisms consisting of two cell types a 1 and a 2 (the case of a greater number of cell types can be studied in a similar way). We can encode these binary cell types setting a 1 = 0 and a 2 = 1. We follow classical ideas [10,9]. The cell pattern is a result of terminal differentiation, which goes by morphogens. Suppose w(x, y, t) is the concentration of a morphogen. Following the positional information concept, one can assume that the cell type depends on the morphogen level. We can suppose that the cell-centered atx j accepts the state a l at the moment t k if open sets in R, which are subsets of the set of possible values of w.
Let t 1 < t 2 , ... < t K be cell differentiation moments. For simplicity and taking into account the periodicity of internal cell dynamics, we assume that t k = k∆t, k = 1, ..., K. It is convenient to introduce the following auxiliary functions of t: Then, to generate a prescribed cell developmental program encoded by the matrix A with binary entries a(l, m) ∈ {0, 1}, we should satisfy the following conditions: where a(l, m) ∈ {0, 1}. We proceed as follows. First, we set where β > 0 and w satisfies boundary conditions (A5), (A6). Taking into account asymptotics (A16) for u we obtain (up to small terms) that for large t where C γ,ϵ is a constant. Using the last relation let us compute Z j (t) defined by (A28). We obtain then Suppose β >> 1 and β << ϵ −1/2 . The points x j,k , which lie in the same k-th group G k , are ϵ 1/2 -close each to the other. On the other hand, the distances between different groups G k have the order O(1) and do not depend on ϵ. We choosex j such thatx j is close enough to j-th group G j and dist(x j , G k ) = O(1) for j ̸ = k. Then for small ϵ > 0 we can simplify (A32) that gives Let us fix j. The dynamics of X l,j (t) from j-th group is defined by j-th system (A27). This means that for large times t all the variables X l,j (t) depend on t via q thus Z j (t) ≈ Z j (q(t)). Then (A29) can be rewritten as where and V j,m are some subsets of the ball B 3 R0 having non-zero Lebesgue measures. These sets are the inverse images (preimages) of the sets W a(j,m) under the maps q → Z j (q).
In the next subsection, we will how to satisfy this relation (A34) and (A35) by chaotic dynamics.

Cell pattern algorithm
An algorithm to solve the pattern generation problems is based on celebrated results of dynamical system theory on hyperbolic sets, its persistence, the existence of Markov partitions and a connection between Bernoulli shifts and maps on invariant hyperbolic sets [11,12,13] can be found in [14].
Let us consider an important particular case of the Markov dynamics described above, when the Markov chain (process) is defined by the system of stochastic differential equations where Q, B are smooth vector field defined on a ball B n in R n , ξ = (ξ 1 , ..., ξ p ) are standard independent white noises. We suppose that the n-dimensional field Q is directed toward to the interior of B n at its boundary, and B is a pdimensional field vanishing at the boundary of that ball. This chain is a process with continuous time, but we can transform it into the discrete-time chain by the Poincaré sections transversal to the flow of the system (A36). Then we obtain a randomized Poincaré map. We suppose that for B = 0 (in the absence of noise) eq. (A36) defines a smooth dynamical system defined with an attractor Γ , which has an invariant measure µ defined on Γ (for axiom A attractors and Anosov diffeomorphisms such measures exist and they are well studied, they are called Sinai-Ruelle-Bowen (SRB) measures [15,16]). Suppose that the flow S t generated by system (A36) has a strong mixing property, i.e.
for two measurable sets A, B. Let µ(A 0 ) be a subset of non-zero µ-measure on Γ and V 0 be a small open neighborhood of A 0 in B n . Let E 1 , E 2 , ..., E nc be a fixed partition of Γ and p be a fixed positive integer. Then the mixing property implies that if ∆T is large enough all the following intersections are non-empty: In fact, according to the strong mixing property µ(B kj ) > 0. Moreover, it is easy to see that This shows that for any finite sequence {a j }, j = 1, ..., p of a j ∈ K c there is an open (possibly having a small measure) set of initial points X(0) such that S j∆T (X(0)) ∈ E aj , j = 1, ..., p.
Let us note that the sophisticated algorithm using Poincaré maps proposed in [14] has an advantage with respect to the simplified one described here: by the Bernoulli shifts and the Markov partitions, we can find the set of initial data q(0) in an explicit way. For our simplified algorithm it can be done by numerical simulations. We have checked it for the Lorenz system for q.
As an example of the simplified algorithm application, let us consider first how to create the French flag pattern (the standard illustration of Wolpert gradient approach).
Let all layers of the cell pattern have the same width as the string corresponding to the French flag, which is s = (1, 2, 3). Let us take the noisy perturbed Lorenz system for q = (q 1 , q 2 , q 3 ), which has the form where η i (t) are independent standard Wiener processes. We suppose that ϵ > 0 is a small parameter, which defines the noise level. We make the standard choice of the parameters σ, ρ, b leading to a chaotic Lorenz attractor Γ ( σ = 10, β = 8/3 and ρ = 28 ).
Let us introduce q min = min q∈Γ q 1 , q max = max q∈Γ , and ∆q = q max − q min . Then we make the partition R 1 = {q : q 1 ∈ (q min , (q max − q min )/2)}, R 2 = {q : q 1 ∈ (q min + (q max − q min )/2, q max )}. The cell fate is determined by the averaged signalq at the moments ∆T, 2∆T... when the wave arrives at the cell. This timeaveraging helps to repress fluctuation effects. Our numerical algorithm is as follows. We use the standard Euler method with a small step, where at each step we generate random numbers dη j . Numerical simulations show that there are starting points q(0) and ∆T such that for exampleq 1 (∆T ) ∈ R 1 , q(2∆T ) ∈ R 2 , andq 1 (3∆T ) ∈ R 1 . So, we obtain the layered 1D-pattern consisting of 3 layers encoded by the string 121. The same construction allows us to obtain more complicated patterns, for example, consisting of 5 and more layers. Let us note that one can take other sets R k , however, these sets should intersect the Lorenz attractor so that the RBS measures of intersections are not too small. Note that the longer the cell pattern s becomes, the smaller the set A 0 of appropriate starting points q(0) will be.

Realisation of vector fields (RVF)
The main technical tool in proving attractor complexity for partial differential equations and systems is the realization vector field (RVF) method using, in particular, structural stability ideas. It is based on a classical center manifold technique and on the well known idea that any n-dimensional dynamics can bifurcate from an equilibrium with n-zero eigenvalues if the number of bifurcation parameters is large enough. Such approach was used for finite dimensional systems (see [17]), but it can be extended on infinite dimensional evolution equations. This RVF approach is developed by first P. Poláčik to prove the existence of non-trivial large-time behavior for quasilinear parabolic equations (see [18,19]) and developed in [20,1] for reaction-diffusion systems, in [21] for neural networks and in [22] for weakly compressible Navier-Stokes equations.
In our model, n = N , where N is the number of kinks, eigenvalues are exponentially close to zero and the initial data z 0 plays the role of the main bifurcation parameter, i.e., the bifurcation parameter is infinite-dimensional.

Structural stability
Recall the basic concept of structural stability introduced by A. Andronov and S. Pontryagin in 1937. Consider a smooth vector field Q on a compact domain D n of R n with a smooth boundary (or on a compact smooth manifold M of dimension n). Assume that Q ∈ C 1 (D n ) and consider all δ-small perturbations Q such that Consider systems of differential equations dq/dt = Q(q) and dq/dt = Q(q) + Q(q) and suppose that they define global semiflows S t Q and S t Q+Q on D n . The system dq/dt = Q(q) is called structurally stable if there exists a δ 0 such that if then trajectories of semiflows S t Q and S t Q+Q are orbitally topologically conjugate (there exists a homeomorphism, which maps trajectories of the first system into trajectories of the second one). Roughly speaking, the original system is structurally stable if any sufficiently small C 1 perturbations of that system conserve the topological structure of its trajectories, for example, the equilibrium point stays an equilibrium (maybe, slightly shifted with respect to the equilibrium of a non-perturbed system), the perturbed cycle is again a cycle (maybe, slightly deformed and shifted).
Note that structurally stable dynamics may be, in a sense, "chaotic". There is a rather wide variation in different definitions of "chaos". We restrict ourselves to hyperbolic chaotic sets. Chaotic (no periodic and no rest point) hyperbolic sets occur in some model systems [23,4,11].

RVF for evolution problems in Banach spaces
Let us consider a family of local semiflows S t P in a fixed Banach space B. Assume these semiflows depend on a parameter P ∈ B 1 , where B 1 is another Banach space. Denote by B n (R) the ball {q : |q| ≤ R} in R n , where q = (q 1 , q 2 , ..., q n ) and |q| 2 = q 2 1 + ... + q 2 n . For R = 1 we will omit the radius R, B n = B n (1). Remind that a set M is said to be locally invariant in an open set W ⊂ B under a semiflow S t in B if M is a subset of W and each trajectories of S t leaving M simultaneously leaves W . In this paper, all W are tubular neighborhoods of the balls B n (R), which have small widths. Consider a system of differential equations defined on the ball B n : where Assume the vector field Q is directed strictly inward at the boundary ∂B n = {q : |q| = 1}: Then system (A39) defines a global semiflow on B n . Let δ be a positive number. Definition (realization of vector fields) We say that the family of local semiflows S t P realizes the vector field Q (dynamics (A39)) with accuracy δ > 0 (briefly, δ -realizes), if there exists a parameter P = P(Q, δ, n) ∈ B 1 such that (i) semiflow S t P has a locally invariant in a open domain W ⊂ B and locally attracting manifold M n ⊂ B diffeomorphic to the unit ball B n ; (ii) this manifold is embedded into B by a map where r > 0; (iii) the restriction of the semiflow S t P to M n is defined by the system of differential equations Definition 2. Let Φ be a set of vector fields Q, where each Q is defined on a ball B n , positive integers n may be different. We say that the family F of local semiflows S t P realizes the family Φ if for each δ > 0 and each Q ∈ Φ the field Q can be δ -realized by the family F.
We say that the family of global semiflows F has the property of universal dynamical approximation if that family realizes the set of all C 1 -smooth finite dimensional fields defined on all unit balls B n .
Many systems enjoy the property of universal dynamical approximation, for example, the Lotka-Volterra system with many species, the Hopfield system, a large class of reaction-diffusion systems, and others. In the coming subsection, we will consider an example of applying the RVF method.

Realization of Rössler system by two component reaction diffusion systems
The well-known Rössler system is one of the simplest systems of ODEs exhibiting chaotic large-time behavior. It reads: where a, b, c are parameters. Under a choice of a, b, c, for example, a = 0.2, b = 0.2 and c = 5.7, the solutions of the Cauchy problem for this system are chaotic (it can be checked numerically). Of course, there exist a number of other systems of ODEs generating chaos, for example, the famous Lorenz system. However, in opposite to the Lorenz system, the Rössler one appears in biochemistry. We will see that there exists an algorithm permitting to realize a system of eqs. (A45), (A46), and (A47) via a class of two-component reaction-diffusion systems. A similar algorithm also realizes all quadratic systems of ODEs (for example, the multispecies Lotka-Volterra systems) by these two-component reaction-diffusion systems.
Consider the following equations for two reagents u, v: where the smooth functions W, g, g k , h and ϵ, a 0 , D are parameters P. We consider equations (A48) and (A49) in the strip Π = (−∞, ∞) × [0, L 2 ] under the following boundary conditions Such a choice of boundary conditions simplifies the problem analysis permitting for quasi-one-dimensional layered solutions for u.
Below we outline the proof of this assertion illustrated by some numerical simulations for a finite difference approximation of the system (A48), (A49). The proof follows the ideas stated above and in [8,24,25,1].

Construction of W , g and asymptotic solutions
The first key step is a choice of W (x), where we use ideas well known in the theory of the Schrödinger operators (see, for example, [25]). Let us consider the 1D Schrödinger operator H = ϵ 2 d 2 dx 2 + W (x) on (−∞, ∞). We can adjust such potential W that this operator will have three eigenvalues λ k and three well-localized eigenfunctions ψ k = ψ(ϵ −1 (x −x k )) with 0 <x 1 <x 2 <x 3 < L 1 corresponding to a discrete spectrum such that λ k = 0, k = 1, 2, 3 whereas all the remaining spectrum is separated by a barrier σ > 0 from the imaginary axis and lies in the (−∞, −1), where σ > 0 is independent of ϵ. Note that since we consider the Neumann conditions at y = 0, L 2 , the same eigenfunctions are "almost" eigenfunctions of the two-dimensional operator ϵ 2 ∆ + W (x), which satisfies periodic boundary conditions with an exponentially small accuracy O exp(−cϵ −1 ) . It permits us to use them to construct an asymptotic solution (for technical details, see [1]).
The solution u can be then represented as where X i are slowly evolving amplitudes andũ is a small correction orthogonal to ψ k in L 2 (Π). The standard ideas for slow-fast decomposition show then [8,24,1] whereṽ is a small correction, and the main contribution V can be represented as sum of three terms: Let us assume that ||ψ k || = 1 and they are orthonormal. The time evolution of X k is governed by a shorted Galerkin system: where (f, g) denotes the standard inner scalar product in L 2 (Π). We can simplify (A55) by substituting asymptotics (A52) and (A53) into the right-hand side of (A55). Taking into account that ψ k are localized in ϵ-neighborhoods at x k we obtain (using the Fourier decomposition) the following asymptotics for V k : whereĝ mk are the Fourier coefficients of g(x k , y): Then we obtain that (A55) reduces to the following quadratic system: where k = 1, 2, 3, functions r k (X, ϵ) are small corrections and where c 0 , c 1 are positive constants depending on the function ψ(z). We obtain the Rössler system if (A61) Therefore, the problem of realization reduces to the following: to adjust the Fourier coefficientsĝ k of g, and the functions h, g 1 , g 2 so that the relations (A60) and (A61) hold. Due to the particular structure of these relations, this problem can be solved in a simple way. Let the support of h, g 2 be disjoint with intervals (−δ +x 1 , δ +x 1 ) and (−δ +x 2 , δ +x 2 ) and the support of a smooth g 1 is disjoint with (−δ +x 3 , δ +x 3 ). This choice provides that h 1 = 0, h 2 = 0 and B ij = 0 for all j = 1, 2, i = 1, 2. Then relations (A60) and (A61) reduce to the following independent linear algebraic systems with respect to unknown Fourier coefficients: where C ij , i = 1, 2, 3, j = 1, 2, 3 are given numbers and b m = ((−1) m − 1)/(mπ/L 2 ) are integrals of sin(myπ/L 2 ) over [0, L 2 ], which do not equal zero for odd integers m. For each j the corresponding system consists of three equations with a countable set of unknownsĝ mj . It is natural to seek solutions with the minimal norm m |ĝ mj | 2 . The numerical solution (see the next subsection) shows that such solutions can be obtained as follows. For each j we choose three odd integers m, say, m 1j = 1, m 2j = 3, m 3j = 5 or m 1j = 1, m 2j = 3, m 3j = 7. We setĝ m,j = 0 for all m ̸ = m lj and we are seeking for the coefficients y l,j =ĝ m l j by solving the linear algebraic systems with respect to y l,j =ĝ m lj j , l = 1, 2, 3. Comments.
(a) For the Rössler system, the construction is simpler than in the general case, although based on the same ideas. For example, the Lorenz system is much more difficult to realize via reaction-diffusion systems. It is connected with the fact that the Lorenz system has no biochemical interpretation, in opposition to the Rössler system.
(b) The chosen boundary conditions simplify the mathematical analysis because the Neumann conditions for u allow us to construct layered quasi-onedimensional solutions for u. On the other hand, for v it is simpler to use the zero Dirichlet conditions at y = 0, L 2 . In fact, for the zero Neumann conditions the corresponding eigenfunctions have the form cos(mπy/L 2 ) and the averages of these functions along [0, L 2 ] equal zero. Then we are not capable to obtain arbitrary matrices A, B.
(c) The difference between this model and the classical approach starting with seminal Turing paper [10] is as follows. The classical morphogenesis theory for two-component models is based on the concept of the Turing pair, one reagent is an activator and the other is an inhibitor.
In this model, which involves spatial inhomogeneities, one of the reactants is neither an activator nor an inhibitor (it is a necessary condition to obtain a stable chaotic behavior in two-component systems, see [1] for mathematical details). To understand this situation, let us consider, for example, eq. (A49) and the reaction term g(x, y)u, which defines the production of v. The function g has a complicated form and changes the sign. In space domains, where g > 0, the reagent u serves as an activator, in domains, where g < 0, u works as an inhibitor for v. Similarly, for the equation (A48) with the term (W (x) + g 2 v)u. Here the reagent u is an activator in layers where W (x) > −g 2 v and it is an inhibitor in other layers. So, we have a mosaic pattern of inhibition and activation. Such structure is necessary to provide a complex non-local interaction between X 1 , X 2 and X 3 , which can produce chaos (see the next comment).
(d) Physical principles are beyond mathematics. The physical and biological ideas are quite standard and follow [26], where is shown a key role of local self-enhancement and long-range inhibition. In our model, local self-enhancement generates narrow vertical layers, where u is concentrated. In our model, we have three layers. The average u concentrations in layers are proportional to X 1 (t), X 2 (t) and X 3 (t). The layers are very narrow thus we have no essential direct interaction between layers. However, we have a long-range reagent v. The u acts to v creating the smooth pattern V (x, y, X 1 , X 2 , X 3 ) (see Figs A1, A2). The new idea with respect to [26] is only the following. The mosaic structure of inhibition and activation induced by spatial gradients (described above) allows us to create complicated dynamics of X i .
(e) Why there is an exponential number of local attractors? It is simpler to explain by the following simplified model, which is more tractable for numerical simulations and considered in the next subsection.  Making numerical calculations for singularly perturbed systems like (A48), (A49) is not an easy problem [25]. However, we can simplify the problem by replacing the equation (A48) with its simplest approximation on a twodimensional lattice. For such an approximation, instead of sophisticated arguments from the Schrödinger operator theory, we can use a much more simple construction. The equation (A49) for v can be replaced by an integral representation by the Green function. Then we obtain a system of ordinary differential equations on the rectangular grid: (A63) where i ∈ {1, 2, ..., L 1 }, j = {1, 2, ..., L 2 }, and 1, t). From a biological point of view, we can consider these equations as a description of colonies of cells located at fixed positions interacting via a reagent v, which can diffuse between cells. Each cell produces a reagent u, which does not diffuse between layers of cells in x-direction.
For V we use the equation    where G(x, y|x,ỹ) is the Green function for the operator D∆ − a 2 0 under boundary conditions (A51). The choice of f i is simple now. Let 1 < i 1 < i 2 < i 3 < L 1 be three integers, which will be positions of layers where the reagent u is concentrated. We set wheref , λ k are positive and λ k <<f . For h i we set For large times t the solution u has the form In fact, all u(i, j, t) with i ̸ = i k are small due to the choice of f i . The quantities u(i k , j, t) are weakly depend on j if d is large enough. Then for X k we have a system of differential equations with quadratic nonlinearities like to the Rössler system (under an appropriate choice of g, λ and b). Comment. The following parameters play a key role in estimating the attractor number: the length of the system L, diffusion rate D for a fast component, ϵ 2 diffusion rate for a slow component, and a 2 0 -degradation rate for a slow component. The distance r K = D 1/2 /a 0 is a measure of the distance traveled by a molecule during its lifetime. It is sometimes referred to as the For each trial, we construct a random stochastic matrix with an absorbing state and choose a random starting state. Then we find τ and after that, we take the average over all trials.
Kuramoto length. The number of possible attractors can be roughly estimated as exp(L/r K ) if ϵ 2 << D. How to explain it?
The complicated dynamics are generated by non-local interaction of closely located u-layers via fast diffusing component (it is sufficient to have 3 layers according to an example from the newly added subsection). The groups of layers should be separated to provide their almost independent functioning. If for dynamics in each group, we have two different local attractors, then if we have M groups thus the full number of attractors is 2 M . So, we have exponential dependence. The groups are well separated if distances between groups are more than the Kuramoto length. It leads to the formula exp(L/r K ). This number is an estimate from below but it is hard to believe that it may be improved. Note that exponential estimates of attractor numbers can be obtained for other models, such as the NK model.
In real conditions, the diffusion rate is bounded. Therefore, to provide adaptivity, the system should increase its size L. We see that the logarithm of τ is approximately N/2 that corresponds to the average Hamming distance between two random bit stings of the length N . The low line describes the case when we decompose the search into three stages. Both curves are obtained by Monte-Carlo simulations over Ntr = 50 trials. For each trial, we proceed random walk on a hypercube choosing a random initial state and a random target state. For each trial we find τ and after we take the average over all trials.
1.10 Some numerical results for the system (4) We denote W the n × n transition matrix with the entries w(s, s ′ ) and by p the vector of state probabilities. Then iterations (4) satisfy the matrix equation p(t + ∆t) = W(ξ(t))p(t). If ξ is absent and the transition probabilities are constants then the mean absorption time τ s to reach the absorbing state from initial state s can be found as a positive solution of the linear algebraic system In the general case of the entries, w depending on ξ we can obtain estimate τ s >τ s , whereτ s are absorption times for the Markov chain with the constant transition probabilities defined by the relationsw(s, s ′ ) = min ξ w(s, s ′ , ξ) for s, s ′ ̸ = s * andw(s * , s) = 1 − s ′w(s * , s ′ ).
We simulate this Markov evolution process by taking W with random entries. Results are represented in Figs. A7 and A9.